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Abstract. We present a new flexible wavefront propagation algorithm for the boundary value problem for sub- 
Riemannian (SR) geodesics in the roto-translation group SE( 2) = R 2 x S 1 with a metric tensor 
depending on a smooth external cost C : SE( 2) —>• [5, 1], S > 0, computed from image data. The 
method consists of a first step where a SR-distance map is computed as a viscosity solution of 
a Hamilton-Jacobi-Bellman (HJB) system derived via Pontryagin’s Maximum Principle (PMP). 
Subsequent backward integration, again relying on PMP, gives the SR-geodesics. For C = 1 we 
show that our method produces the global minimizers. Comparison with exact solutions shows a 
remarkable accuracy of the SR-spheres and the SR-geodesics. We present numerical computations 
of Maxwell points and cusp points, which we again verify for the uniform cost case C = 1. Regarding 
image analysis applications, tracking of elongated structures in retinal and synthetic images show 
that our line tracking generically deals with crossings. We show the benefits of including the sub- 
Riemannian geometry. 

Key words. Roto-translation group, Hamilton-Jacobi equations, Vessel tracking, Sub-Riemannian geometry, 
Morphological scale spaces. 
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1. Introduction. In computer vision, it is common to extract salient curves in images 
via minimal paths or geodesics minimizing a length functional. The minimizing geodesic is 
defined as the curve that minimizes the length functional, which is typically weighted by a cost 
function with high values on image locations with high curve saliency. To compute such data- 
driven geodesics many authors use a two step approach in which firstly a geodesic distance 
map to a source is computed and then steepest descent on the map gives the geodesics. In 
a PDE framework, the geodesic map is obtained via wavefront propagation as the viscosity 
solution of a Hamilton-Jacobi-Bellman equation (the Eikonal equation). For a review of this 
approach and applications see [37, 27, 31]. 

Another set of geodesic methods, partially inspired by the psychology of vision was de¬ 
veloped in [15, 30]. Here, the roto-translation group SE( 2) = M 2 x S 1 endowed with a 
sub-Riemannian (SR) metric models the functional architecture of the primary visual cortex 
and geodesics are used for completion of occluded contours. A stable wavelet-like approach 
to lift 2D-images to functions on SE( 2) was proposed in [16]. Within the SE( 2) framework, 
images and curves are lifted to the 3D space M 2 x S 1 of coupled positions and orientations in 
which intersecting curves are disentangled. The SR-structure applies a restriction to so-called 
horizontal curves which are the curves naturally lifted from the plane (see Fig. 1A). For ex¬ 
plicit formulas of SR-geodesics and optimal synthesis see [36]. SR-geodesics in SE( 2) were also 
studied in [8, 10, 18, 23, 25]. Here, we propose a new wavefront propagation-based method 
for finding SR-geodesics within SE( 2) with a metric tensor depending on a smooth external 
cost C : SE( 2) -a [J, 1], S > 0 fixed. Our solution is based on a Hamilton-Jacobi-Bellman 
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Figure 1 . A: Every point in the planar curve 72D (t) = (x{t),y(t)) is lifted to a point g = 7 (t) = 
(x{t),y(t), 0 (t)) G SE( 2 ) on an horizontal curve (solid line) by considering the direction of the tangent vector 
72D (t) of the planar curve as the third coordinate. Then, tangent vectors 7 (t) G span{ Ai | 7 ^ , A2\ 1 ^} = A| 7 ^ ; 
Eq. ( 2 . 1 ). B: In the lifted domain SE( 2) crossing structures are disentangled. C: The SR-geodesic (green) better 
follows the curvilinear structure along the gap than the Riemannian geodesic (red). 


(HJB) equation in SE( 2) with a SR metric including the cost. It is of interest to interpret 
the viscosity solution of the corresponding HJB equation as a sub-Riemannian distance [38]. 
Using Pontryagin’s Maximum Principle (PMP), we derive the HJB-system with an Eikonal 
equation providing the propagation of geodesically equidistant surfaces departing from the 
origin. We prove this in Thm. 3.1, and we show that SR-geodesics are computed by back¬ 
tracking via PMP. In Thm. 3.2, we consider the uniform cost case (i.e. C — 1) and we show 
that the surfaces coincide with the SR-spheres, i.e. the surfaces from which every tracked 
curve is globally optimal. This uniform cost case has been deeply studied in [36] relying on 
explicit ODE-integration in PMP. In this article, we will rely on a PDE-approach and viscosity 
solutions of HJB-equations, allowing us to extend to the general case where C is a smooth cost 
uniformly bounded from below and above. We will often use the results in [36] as a golden 
standard to verify optimality properties of the viscosity solutions and accuracy of the involved 
numerics of our PDE-approach. We find a remarkable accuracy compared to exact solutions, 
1st Maxwell sets (i.e. the location where for the first time two distinct geodesics with equal 
length meet), and to the cusp surface [18, 10]. 

Potential towards applications of the method with non-uniform cost is demonstrated by 
performing vessel tracking in retinal images. Here the cost function is computed by lifting 
the images via oriented wavelets, as is explained in Section 4.1. Similar ideas of computing 
geodesics via wavefront propagation in the extended image domain of positions and orienta¬ 
tions, and/or scales, have been proposed in [28, 21, 7]. In addition to these interesting works 
we propose to rely on a SR geometry. Let us illustrate some key features of our method. In 
Fig. IB one can see how disentanglement of intersecting structures, due to their difference in 
orientations, allows to automatically deal with crossings (a similar result can be obtained with 
the algorithm in [28]). The additional benefit of using a SR geometry is shown in Fig. 1C 
where the SR-geodesic better follows the curvilinear structure along the gap. 

1.1. Structure of the Article. The article is structured as follows. First, in Section 2, we 
give the mathematical formulation of the curve optimization problem that we aim to solve in 
this paper. We also provide embedding into previously studied geometrical control problems 
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formulated on SE( 2) and/or R 2 . In Section 3 we describe our PDE approach that provides 
the SR geodesic distance map as the viscosity solution of a boundary value problem (BVP) 
involving a sub-Riemannian Eikonal equation. Furthermore, in Theorem 3.1, we show that 
sub-Riemannian geodesics are obtained from this distance map by back-tracing (imposed by 
the PMP computations presented in Appendix A). In Theorem 3.2 we show that for the 
uniform cost case (i.e. C — 1) such back-tracing will never pass a Maxwell-points or conjugate 
points, and thereby for C = 1 our approach provides only the globally optimal solutions. 

In Section 4 we describe an iterative procedure on how to solve the BVP by solving a 
sequence of initial value problems (IVP) for the corresponding HJB-equation. Before involve¬ 
ment of numerics, we express the exact solutions in concatenated morphological convolutions 
(erosions) and time-shifts in Appendix E. Here we rely on morphological scale space PDE’s 
[12, 2, 17], and we show that solutions of the iterative procedure converges towards the sub- 
Riemannian distance map. Then in Section 4.1 we construct the external cost C, based on a 
lifting of the original image to an orientation score [16]. 

In Section 5, we describe a numerical PDE-implementation of our method by using left- 
invariant finite differences [20] in combination with an upwind-scheme [33]. Finally, in Sec¬ 
tion 6, we present a verification of the proposed method by comparison with the exact solu¬ 
tions. We also provide simple numerical approaches (extendable to the non-uniform cost case) 
to compute 1st Maxwell points, conjugate points and cusp points [18], which we verify for the 
uniform cost case with results in [36]. Finally, in Subsection 6.2, we show application of the 
method to vessel tracking in optical images of the retina. We discuss the two main parameters 
that are involved: the balance between external and internal cost, and the balance between 
spatial and angular motion. First feasibility studies are presented on patches, and we discuss 
on how to proceed towards automated retinal vessel tree segmentation. 

2. Problem Formulation. The roto-translation group SE( 2) carries group product: 
gg' = (x, Re)(x.', Re') = (Rex' + x, R e+ e')- 

where Rq is a counter-clockwise planar rotation over angle 9 . This group can be naturally 
identified with the coupled space of positions and orientations R 2 x S' 1 , by identifying Rq ^ 9 
while imposing 27r-periodicity on 6. Then for each g G SE( 2) we have the left multiplication 
L g h = gh. Via the push-forward ( L g )* of the left-multiplication we get the left-invariant 
vector fields {Mi, M2, M3} from the Lie-algebra basis {Mi, M2, M3} = {d x \ e , do\ e , d y \ e } at the 
unity e = (0,0, 0): 


Ai\ g = cos(9 d x \ g + sin<9 d y \ g = (L g )* d x \ e , 

(2M) M 2 |^ = do\ g = (L g )* do\ e , 

M 3 | g = - sin 9 d x \ g + cos 9 d y \ g = (L g )* d y \ e . 

Then all tangents y(t) G T^(SE( 2)) along smooth curves t y(t) = (#(£), y(£), 0(£)) G 
SE( 2) can be expanded as y(t) = Ylk=i uk (f) M&| 7 (fp where the contravariant components 
u k (t ) of the tangents (velocities) can be considered as the control variables. Not all curves t 1 — 
7 (t) in SE( 2) are naturally lifted from the plane in the sense that 9{t) = aig(x(t)+i This 

holds for so-called horizontal curves which have u 3 = 0 and thus y(t) = Y^k=i uk (t) M&| 7 (t)- 
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The allowed (horizontal) directions in tangent bundle T(SE( 2 )) (see Fig. 1 A) form a so-called 
distribution 

A := span{*4i,*4 2 }- 

Therefore we consider SR-manifold [26] (SE{2), A, G c ), with 
G c : SE( 2 ) xAxA^R denoting the inner product given by 

( 2 . 2 ) G c | 7 (t)( 7 (t), 7 (i)) = C 2 (j(t)) (p 2 \x(t)cos0(t)+y(t)sm6(t)\ 2 + |0(t)| 2 ) , 

with 7 : 'R—>SE( 2 ) a smooth curve on R 2 xi S' 1 , (5 > 0 constant, C : SE( 2 ) —X [5,1] the given 
external smooth cost which is bounded from below by 6 > 0. 

Remark 1. Define C g (j){h ) = </>(g - 1 /i) t/ien iee have: 

= G c ° c \ gi UL g )^,(L g )^). 

Thus, G c is not left-invariant, but if shifting the cost as well, we can, for the computation of 
SR-geodesics, restrict ourselves to 7 ( 0 ) = e. 

We study the problem of finding SR minimizers, i.e. for given boundary conditions 7 ( 0 ) = 
e, 7 (T) = g , we aim to find the horizontal curve 7 (t) (having 7 E A) that minimizes the total 
SR length 

(2.3) 1=1 y/G c | 7(t) Cf(t),j(t))dt. 

a 0 

If t is the SR arclength parameter, our default parameter, then G c |, (7(f), 7 (t)) = 1 and 

l = T. Then, SR minimizers 7 are solutions to the optimal control problem (with free T > 0): 


(2.4) 


PLc(S£(2)) : ^ 


7 = u l Ai| 7 + u 2 A 2 \ 7 , 

7(0) = e, 7CO = 9, 

K'y(')) = lo c {lX))\/WWT)? + |« 2 (i)| 2 di —> min, 
7 (f) € SE( 2 ), (u 1 ^), u 2 (t)) € R 2 , f3 > 0. 


In the naming of this geometric control problem we adhere to terminology in previous work 
[10, 18]. Stationary curves of the problem (2.4) are found via PMP [1]. Existence of minimizers 
follows from Chow-Rashevsky and Filippov’s theorem [1], and because of absence of abnormal 
trajectories (due to the 2 -bracket generating distribution A) they are smooth. 

Remark 2. The Cauchy-Schwarz inequality implies that the minimization problem for the 
SR length functional l is equivalent (see e.g. [26]) to the minimization problem for the action 
functional with fixed T: 

J (7) = \J Q c 2 Ci(t)){P 2 \ ul X )\ 2 + l« 2 WI 2 ) dt - 


(2.5) 
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2.1. Embedding into Geometric Control Theory. The problem (2.4) actually comes from 
a mechanical problem in geometric control, where a so-called Reeds-Shepp car [32] can proceed 
both forward and backward in the path-optimization. As pointed out in [9] such a problem, for 
certain end conditions, cannot be considered as a curve optimization problem on the plane. In 
fact in [17] the set 94 C SE( 2) of allowable end-conditions g is explicitly determined for C — 1. 
The boundary of the set 94 is formed by end-points of geodesics whose spatial projections 
either depart or end in a cusp [17, Fig.l2,Thm.9]. For g — (x, 6) G 94 C SE( 2) and C — 1 the 
following problem is well-posed 


' 7 ( 0 ) = 0, 7 (L) = x, 

7(0) = (r°) T , 7 ( L ) = (cos 9, s'mO) T , 

L 

l(j(-)) = J ^(l( s )) + k 2 (s) d s -> min, 

. 7 : [0, L] —> R 2 , (5 > 0, 


where L denotes spatial length and k curvature of smooth curve 7 € G'°°([ 0 , L] , M 2 ). For C / 1 
the set of end-conditions (x, 9) for which the above problem is well-posed varies with C. 

In this article we take P^ nec (SE(2)) (and not (2.6)) as a venture point, as for the detection 
of elongated structures we must look both forward and backward. From the image analysis 
point of view, our primary interest is in the following refined problem of P^ iec (S f F^(2)): 


(2.7) 


Pco„tour(S£(2)) : < 


7 (<) = ^(t) + u 2 {t) v4 2 | 7 p) , for t G [0, T] 

7 (0) =e, 7(0 = g = (x,0) € Dl c , 

l(j(-)) = J C( 7 (t))-v// 3 2 |n 1 (t )| 2 + |?/ 2 (t)| 2 dt — >■ min, 
0 

with curve 7 : [0,T] SE(2), with controls: 

(u 1 (t), u 2 (t)) G M 2 , and u l {t) does not change sign, 


where 94 c is the set of all g G SE( 2) such that the minimizing SR-geodesic(s) y(-) = (x(-), O(-)) 
do not exhibit a cusp in their spatial projections x(-). 

Remark 3 .For existence results of the minimizers, a priori Li -conditions have to be imposed 
on the controls u\, U 2 . However, in retrospect, after application of a generalized version of 
PMP [39], [10, Thm 5.2 and App.B], the controls of the minimizers turn out to be smooth 
and bounded. Regarding extension to C ^ 1 (recall (2.2)) we note that S < C < 1 . 

Remark 4 .For C — 1 one has 


7f c=1 =9luO with Q = {(x, y, 6) € SE( 2) | (-x, y, -0) € 91}, 

and with 94 C {(x,y,6) G SE( 2) | x > 0} described explicitly in [18], and partly depicted in 
item C of Figure 2. In Section 6.1 we shall provide a very simple numerical tool to compute 
the surface in SE( 2 ) where cusps appear also for C ^ 1 . This surface is a boundary of a 
volume in SE( 2 ) that contains the set 94 c . 

Summarizing, in P ( [ nec (SE(2)) the optimal control t u x {t) can switch sign. 

• If g is chosen such that the optimal control u 1 > 0 then the lift of problem Pcurve(^ 2 ) 
coincides with P^ iec (S'£’(2)) and also with Pcontour(^^(2))- 
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Figure 2. 3 Sub-Riemannian geodesics for uniform cost C = 1 to reveal the differences of the 3 geometric 
control problems (2.4), (2.6) and (2.7), for /3 = 1. A: plots of spatial projections of sub-Riemannian geodesics 
in R 2 , B: SR-geodesics in SE( 2), C: A part of (the set of end-conditions for which P CU rve{ R 2 ) is well- 
defined [18]) depicted as reachable cones around the origin. End condition g\ — (x\,yi,6i) G yields the 
minimizing curve of P contour E 2)); Pmec (SE( 2)) and P CU rve( R 2 ). End condition g 2 = (— an, yi, — 0i) yields 
the minimizing curve in P con tour{SE{2)), P mec (SE( 2)), and it is invalid for P curt;e (R 2 ). End condition gz is 
invalid for both P CU rve( R 2 ) and P C ontour(SE(2)), as it induces an internal cusp. For C — 1 the set of allowable 
end conditions for P CO ntour(SE(2)) equals 9t c=1 = 9tU 0,. For C ^ 1 this set 9t c differs, and can be computed. 


• If g is chosen such that the optimal control u 1 < 0 then problem P^ iec (S f £^(2)) and 
problem Pcontour(*S^(2)) coincide. 

• If g is chosen such that the optimal control u 1 (t) switches sign at some internal time 
t G (0,T), then g G SE( 2) \ and the spatial projection of the corresponding 
minimizing SR-geodesic(s) has an internal cusp, which we consider not desirable in 
our applications of interest. 

See Figure 2, where each of the 3 above cases is illustrated. 

Now let us return to our main objective, that is to derive solutions of (2.4) via a PDE- 
wavefront propagation method for data-driven SR-geodesics. The word ‘data-driven’ refers to 
the fact that our PDE-approach is also applicable to C ^ 1, a property that does not apply 
to previous methodology following ODE-approaches [36, 10, 23, 18, 30] on SR-geodesics in 
vision/image analysis. 

3. Solutions via Data-driven Wavefront Propagation. The following theorem summa¬ 
rizes our method for the computation of data-driven sub-Riemannian geodesics in SE( 2). The 
idea is illustrated in Fig. 3. 

Theorem 3.1. Let W(g ) be a solution of the following boundary value problem (BVP) with 
Eikonal- equation 


(3.1) 


V(C0))- 2 + \AWW) -1 = 0, for g^e 

W(e) = 0. 
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Then the iso-contours 

(3.2) S t = {g € SE( 2 ) | W(g) = t} 

are geodesically equidistant with speed = C(^(b))yJ fi 2 \u l {t)\ 2 + \u 2 {t)\ 2 = 1 and they provide 
a specific part of the SR-wavefronts departing from e = (0, 0, 0). A SR-geodesic departing from 
g G SE( 2) is found by backward integration 

( 3 - 3 ) 7fe ( ) (/?C( 7b (i))) 2 (C( 76 (t))) 2 7 ^°) 9 ■ 

Proof The definition of geodesically equidistant surfaces is given in Definition B.l in Ap¬ 
pendix B. Furthermore, in Appendix B we provide two lemmas needed for the proof. In 
Lemma B.2, we connect the Fenchel transform on A, to the Fenchel transform on R 2 to ob¬ 
tain the result on geodesically equidistant surfaces in (SE(2), A, G c ). Then, in Lemma B.3 
in Appendix B, we derive the HJB-equation for the homogeneous Lagrangian as a limit from 
the HJB-equation for the squared Lagrangian. The back-tracking result follows from applica¬ 
tion of PMP to the equivalent action functional formulation (2.5). Akin to the R^-case [11], 
characteristics are found by integrating the ODE’s of the PMP where according to the proof 
of Lemma B .2 we must set p = d SR W, see Remark 5 below. □ 

The next theorem provides our main theoretical result. Recall that Maxwell points are 
SE( 2 ) points where two distinct geodesics with the same length meet. The 1st Maxwell set 
corresponds to the set of Maxwell-points where the distinct geodesics meet for the first time. 
In the subsequent theorem we will consider a specific solution to (3.1), namely the viscosity 
solution as defined in Definition C.3 in Appendix C. 

Theorem 3.2. Let C — 1. Let W(g ) be the viscosity solution of the BVP (3.1). Then St, 
Eq. (3.2), equals the SR-sphere of radius t. Backward integration via (3.3) provides globally 
optimal geodesics reaching e at t = W(g) = d(g , e) := 

(3.4) min [ \f\6{r)\ 2 + f3 2 \x(r) cos 9{r)+y{r) sin 6{r) | 2 dr, 

7 G C°°(M + , 5'E(2)), T > 0 , Jo v 

7 G A, 7 ( 0 ) = e, 7 (T) = g 

and 7 b(t) = 7 mm (d(g,e) — t). The SR-spheres St = {g G SE( 2) | d(g,e) = t} are non-smooth 
at the 1st Maxwell set Ai, cf. [36], contained in 

(3.5) M C { (x, y , 6) G SE( 2) | x cos | + y sin | = 0 V 9 — 7 r} , 

and the back-tracking (3.3) does not pass the 1st Maxwell set. 

Proof of Thm. 3.2 can be found in Appendix D. The global optimality and non-passing 
of the 1st Maxwell set can be observed in Fig. 4. For the geometrical idea of the proof see 
Fig. 5. 

Remark 5. The stationary solutions of (4-1) satisfy the SR-Eikonal equation (3.1). The 
Hamiltonian H^ lxed for the equivalent fixed time problem (2.5) equals 

H fixed (g,p ) = (CO))" 2 {r 2 hl + hi) = 1/2, 


(3.6) 
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Figure 3. A-B: Our method provides both geodesically equidistant surfaces St (3.2) (depicted in A) and 
SR-geodesics. As depicted in B: geodesic equidistance holds with unit speed for all SR-geodesics passing through 
the surface, see Thm 3.1. Via Thm. 3.2 we have that W(g) — d(g,e) and {<St}t>o is the family of SR-spheres 
with radius t depicted in this figure for the uniform cost case. 


Figure 4. A: SR-sphere St for t — 4 obtained by the method in Thm. 3.1 using C — 1 and as initial 
condition via viscosity solutions of the HJB-equation (4-1) implemented according to Section 5. B: The full 
SR-wavefront departing from e via the method of characteristics and formulae in [25] giving rise to interior folds 
(corresponding to multiple valued non-viscosity solutions of the HJB-equation). The Maxwell set M consists 
precisely of the dashed line on xcos | + y sin | = 0 and the red circles at \9\ = tt. The dots are 2 (of the 4) 
conjugate points on St which are limits of 1st Maxwell points (but not Maxwell points themselves). In B we 
see the astroidal structure of the conjugate locus [35, 13]. In A we see that the unique viscosity solutions stop 
at the 1st Maxwell set. Comparison of A and B shows the global optimality and accuracy of our method at A. 
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Figure 5. Maxwell point g* = (—4, 4,7r/2) (in white) on SR-sphere St (in orange) for C — 1. At g* two 
SR-geodesics 71 ^ 72 with equal SR-length t meet (71 (t) = 72 if)). From left to right: A: projection of 71 and 
72 on the plane (x,y), B: 2D-slices (x — x*, y = y*) of level sets of W(g) with distinguished value W(g) = t 
(again in orange). On top we plotted, the Maxwell point, the intersection of surface xcos | + y sin | = 0 (in 
purple, this set contains a part of the 1st Maxwell set) with the 2D-slices. C: The SR-sphere St in SE(2), 
D: section around g * revealing the upward kink due to the viscosity solution. From this kink we see that the 
tracking (3.3) does not cross a 1st Maxwell point as indicated in red, yielding global optimality in Thm. 3.2. 


with momentum covector p — h\uo l + /i 2 <^ 2 + h^uj 3 expressed in dual basis {uj 1 } 3 =1 given by 

(3.7) {uo l ,Aj) — Sj AA uj l = cos#dx + smddy, uj 2 = dd, cv 3 = — smOdx + cosddy. 

The Hamiltonian H^ ree for the free time problem (2.4) minimizing l equals 

(3.8) H free ( g,p) = _i = o. 

For details see Appendix A and B. Eq. (3.1) can be written as Hf ree (g,p) — 0 with momentum 1 

2 

p = d SR W — oA 

2=1 


Remark 6. SR geodesics loose their optimality either at a Maxwell point or at a conjugate 
point (where the integrator of the canonical ODE’s, mapping initial momentum po and time 
t > 0 to end-point 7 (t), is degenerate [1]). Some conjugate points are limits of Maxwell points, 
see Fig. 4, where the 1st astroidal shaped conjugate locus coincides with the void regions (cf. [3, 
fig-1]) after 1st Maxwell set A4. When setting a Maxwell point as initial condition, the initial 
derivative d SR W\^^ is not defined. Here there are 2 horizontal directions with minimal 
slope, taking these directions our algorithm produces the results in Fig. 5 A and Fig. 14 . 


1 Note that the sub-Riemannian gradient V SR W — G 1 dW — Y2l=i Pi 2 AiWAi, with (3\ = /3, /3 2 — 1, by 

definition is the Riesz-representative (being a vector) of this SR-derivative (being a covector). 
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4. An Iterative IVP-procedure to Solve the SR-Eikonal BVP. To obtain an iterative 
implementation to obtain the viscosity solution of the SR-Eikonal BVP given by (3.1), we rely 
on viscosity solutions of initial value problems (IVP). In this approach we put a connection 
between morphological scale spaces [12, 2], and morphological convolutions with morphological 
kernels, on the SR manifold (SE( 2), A, G c ) and the SR Eikonal BVP. 

In order to solve the sub-Riemannian Eikonal BVP (3.1) we resort to subsequent auxiliary 
IVP’s on SE( 2) for each r G [r n , r n + 1 ], with r n — ne at step n G N U {0}, e > 0 fixed: 


(4.1) 


' = 1 - yJ{C{g))~* {p- 2 \AiWZ +1 (g,r)\* + 

< W* +1 (g,r n ) =WZ(g,r n ) fox g^e, 

< ^n+l( e 5 r ^) 


for n — 1 , 2 ,..., and 

(A 2 ^ / = 1 - y(C( 5 ))- 2 + \A 2 W{{g,r)\*), 

l Wl(g,0) =5M(g), 


for n — 0, where 5^ is the morphological delta given by 




0 if g = e 
oo else. 


Let W^+i denote the viscosity solution of (4.1) carrying the following support 

supp(W^ +1 ) = SE( 2) x [r n ,r n+ i], with r n = ne, 

so in (4.1) at the n-th iteration (n > 1) we use, for g ^ e, the end condition W^{g^r n ) of the 
n-th evolution as an initial condition W^+iig, r n ) of the (n + l)-th evolution. Only for g = e 
we set initial condition B / ^ +1 (e,r n ) = 0. Then we define the pointwise limit 


(4.3) W°°(g) := lim ( lim W^ +l (g,ne)) . 

e—^0 \n—>oo J 

Finally, regarding the application of our optimality results, it is important that each IVP- 
solution Wn+i(g,r) is the unique viscosity solution of (4.1), as then via (4.3) the viscosity 
property for the viscosity solutions of the HJB-IVP problem naturally carries over to the 
viscosity property of the viscosity solutions of system (3.1). Thus we obtain W = W°° as the 
unique viscosity solution of the SR-Eikonal BVP. 

Details on the limit (4.3), which takes place in the continuous setting before numeric 
discretization is applied, can be found in Appendix E. In Appendix E we provide solutions 
of (4.1) by a time shift in combination with a morphological convolution 2 with the corre¬ 
sponding morphological kernel, and show why the double limit is necessary. A quick intuitive 
explanation is given in Figure 6, where we see that for e > 0 we obtain stair-casing (due to 
a discrete rounding of the distance/value function) and where in the limit e ^ 0 the solution 
W°°(g) — W(g) — d(g,e) is obtained. 

2 In fact, an ‘erosion’ according to the terminology in morphological scale space theory, see e.g. [12]. 
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Figure 6. Illustration of the pointwise limits in Eq. (4-3). Top: plot of g i—> lim VK^ +1 (g, r n +i) (from 

n—>oo 

left to right, resp. for e — 1, e — 0.5 and e 4 Q) which is piecewise step-function, see Corollary E.2 in 
Appendix E. Along the red axis {(a:, 0,0) | x G M} we have x = d(g,e). Bottom: the corresponding graph of 
x i —y Wn+i((x, 0 , 0 ), r n +i). As n grows the staircase grows, as e —»> 0 the size of the steps in the staircase 
vanishes and we see W°°(g) = d(g,e) in the most right column. 


Remark 7. The choice of our initial condition in Eq. (4-2) comes from the relation between 
linear and morphological scale spaces [2, 12]. Here, for linear SE (2) -convolutions over the 
(-,+)-algebra one has S e *se( 2) U — U. For morphological SE (2) -convolutions (erosions) over 
the (min, +)-algebra [18] one has a similar property: 

(4.4) (<5 e M © 17)0/) := inf {5^{ q ~ l g) + U(q)} = U(g), 

qeSE( 2) 

This is important for representing viscosity solutions of left-invariant HJB-equations on SE( 2) 
by Lax-Oleinik [19] type of formulas (akin to the SE(3)-case [17]). 

Remark 8. The stair-casing limit depicted in Figure 6 is similar to the basic Eikonal BVP 
on M with solution d(x, 0) = \x\. On R the approach (4-1),(4-2) and (4-3) provides pointwise 

oo r^i 

limit |x| = lim ]T r m+1 l [rm , Pra+e] (|*|) = lim £ 1 [r m ,r m + e ](l*D- 

e40 m=0 m=0 

Remark 9 .By general semigroup theory [2], one cannot impose both the initial condition 
and a boundary condition W e (e,r) — 0 at the same time, which forced us to update the initial 
condition (at g — e) in our iteration scheme (4-1)- The separate updating with value 0 for 
g — e in Eq. (4-1) is crucial for the convergence in Eq. (4-3). 


4.1. Construction of the Non-uniform Cost. The cost should have low values on locations 
with high curve saliency, and high values otherwise. Based on image / we define the cost- 
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function 5 < C < 1 via 
(4.5) 


C(x, y, 9) 


1 


1 + A 


(AlUf)+[x,y,9) 

ll(-4^/)+l|oo 


pi 


with (5 = and where Uf : SE( 2 ) —> R is a lift of the image, with || • the sup-norm, and 
(A\Uf) + {g) = { (- sin ^ + cos Od y )*U f (. 


X. 77. 0) if I 


is a vessel detector where we use spatially isotropic Gaussian derivatives [20]. The vessel- 
detector gives responses only if there are convex variations orthogonal to the elongated- 
structures of interest in Uf(x,y,0). 

The lifting is done using real-valued anisotropic wavelets 7 /k 

(4.6) Uf(x.,e)= f -x))/(y)dy. 

J R 2 

In this work we use the real part of so-called cake wavelets [16] to do the lifting. These 
wavelets have the property that they allow stable reconstruction and do not tamper data 
evidence before processing takes place in the SE( 2 ) domain. Other type of 2 D wavelets could 
be used as well. In related work by Pechaud et al. [28] the cost C was obtained via normalized 
cross correlation with a set of templates. 

In Eq. (4.5) two parameters, A and p, are introduced. Parameter A can be used to increase 
the contrast in the cost function. E.g., by choosing A = 0 one creates a uniform cost function 
and by choosing A > 0 data-adaptivity is included. Parameter p > 1 controls the steepness of 
the cost function, and in our experiments it is always set to p — 3. 


/ 




Figure 7. Illustration of the cost function C. Left: retinal image patch f. Middle: corresponding function 
Uf (‘invertible orientation score’) using the real part of a cake-wavelet 'if [16]. Right: the cost function C 
computed via Eq. (f.5) visualized via volume rendering. The orange corresponds to locations where C has a low 
value. 

5. Implementation. To compute the SR geodesics with given boundary conditions we 
first construct the value function W°° in Eq. (3.1), implementing the iterations at Eq. (4.1), 
after which we obtain our geodesic 7 via a gradient descent on W°° from g back to e, recall 
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Thm. 3.1 (and Thm. 3.2). Throughout this section, we keep using the continuous notation 
g G SE( 2) although within all numerical procedures g is sampled on the following (27V + 1) x 
(:27V + 1) x (27V#) equidistant grid: 

(5.1) {(zi,yj,O k ) | Xi = s xy i,Xj = s xy j,Q k = ks e , with i,j = -TV, ...TV, k = -7V# + 1,... TV#}, 

with step-sizes s# = s xy = 2 n+i ? with A} TV# G N. As a default we set TV = 60, = 7, 

TV# = 32. The time-discretization grid is also chosen to be equidistant with time steps Ar = e. 

On this grid we compute an iterative upwind scheme to obtain the viscosity solution W e 
at iteration Eq. (4.1). Here we initialize VE e (-,0) = 4 ( ^ jD (-), with the discrete morphological 
delta, given by S MD (g) = 0 if g = e and 1 if g ^ e, and iterate 


f r + Ar) = r) - Ar H f D ree (g, dW*(g, r)) for g + e 

\ VF e (e, r + Ar) =0, 


with free-time Hamiltonian (see Appendix A, Eq. (A.4)) given by 

H f D ree (g,dW £ (g,r)) = ~ l) , 

until convergence. We set Ar = e in Eq. (4.1). In the numerical upwind scheme, the left- 
invariant derivatives are calculated via 

(• AiW £ (g,r )) 2 = (max {At W e (g, r), -Af W e (g, r), 0}) 2 , 

where Af and denote respectively the forward and backward finite difference approxima¬ 
tions of A{. Note that W e in (5.2) is a first order finite difference approximation of W^ +1 in 
(4.1) at time interval r G [ne, (n + l)e] and we iterate until the subsequent L^-norms differ 
less than 10 -6 . This upwind scheme is a straightforward extension of the scheme proposed in 
[33] for HJB-systems on M n . It produces sharp ridges at the 1st Maxwell set (cf.Fig.4) as it is 
consistent at local maxima. For numerical accuracy and left-invariance we applied finite differ¬ 
ences in the moving frame of left-invariant vector fields, using H-spline interpolation. This is 
favorable over finite differences in the fixed coordinate grid {x, y, 0}. For details on these kind 
of left-invariant finite differences, and comparisons to other finite difference implementations 
(e.g. in fixed coordinate grid) see [20]. 

In our implementation the origin e is treated separately as our initial condition is not 
differentiable. We apply the update W e (e,r) = 0 for all r > 0. We set step size e = 
0.1 min (s xy /3,so) with s xy and sq step sizes in respectively the x-?/-directions and 0-direction. 

6. Experiments and Results. 

6.1. Verification for the Uniform Cost Case. Throughout the paper we have illustrated 
the theory with figures obtained via our new wavefront propagation technique. In this sub¬ 
section we go through the figures that support the accuracy of our method. As the problem 
(2.4) for C = 1 was solved [36, 18], we use this as a golden standard for comparison. 
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Figure 8. A: SR-geodesic example for the uniform cost case shows or PDE-discretizations (with 12 and 64 
sampled orientations in red and green respectively) are accurate in comparison to analytic solutions in [36, 18] 
(in black). B: The blue surface represents the cusp surface numerically computed via the proposed HJB-system 
(with C — 1) and subsequent calculation of the zero-crossings of AiW°°(x,y, 6). Indeed if a SR-geodesic (in 
green) passes this surface, it passes in 6-direction (with infinite curvature [10, 18]), yielding a cusp on the 
spatial ground plane. The same blue surface is computed in [18, Fig. 11]. We even see the additional fold (top 
left passing the grey-plane) as some globally optimal SR-geodesics even exhibit 2 cusps. 


6.1.1. Comparison of BVP Solutions and the Cuspsurface. Let us consider Fig. 8 A. 
Here an arbitrary SR-geodesic between the SE( 2 ) points 7 ( 0 ) = e and 7 (T) = ( 6 , 3, 7 t/ 3 ) is 
found via the IVP in [36] with end-time T — 7.11 and initial momentum 

po = hi(0)dx + /i 2 ( 0 )dp + /i 3 (O)d0, 

with hi(0) — \fl — |^2 (0)| 2 , / 12 (C)) = 0.430 and hs(0) = —0.428, is used for reference (in black 
in Fig. 8 A . Using the semi-analytic approach for solving the BVP in [18] an almost identical 
result is obtained. The curves computed with our method with angular step-sizes of 27r/12 
and 27r/64 are shown in Fig. 8 A in red and green respectively. Already at low resolution we 
observe accurate results. In Fig. 4 we compare one SR-sphere for T — 4 (Fig. 4A) found via 
our method with the SR-wavefront departing from e (Fig. 4B) computed by the method of 
characteristics [25]. We observe that our solution is non-smooth at the 1st Maxwell set M 
(3.5) and that the unique viscosity solution stops precisely there, confirming Theorem 3.2. 
Finally, the blue surface in Fig. 8 B represents the cusp surface, i.e. the surface consisting of 
all cusp points. Cusps are points that can occur on geodesics when they are projected into 
the image plane (see Fig. 8 B). This happens at points g c where the geodesic is tangent to 
do\g c = * 4 . 2 \g c • Then, the cusp surface & CU sp is easily computed as a zero-crossing: 

( 6 . 1 ) 6 ctisp := {g € SE( 2 ) | (A^g, e) = 0 }, 6 ™™ := {g € SE( 2 ) | AiW°°(g) = 0 }. 

It is in agreement with the exact cusp surface & CUS p analytically computed in [18, Fig. 11]. 

Remark 10 .The simple key geometric idea behind (6.1) is that we have a cusp at time t if 
u l (t) = c 2 (^(t)) = = 0 which directly follows from Appendix A. 
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6.1.2. Comparison and Computation of SR-spheres. Numerical verification of the solu¬ 
tion obtained by our PDE approach was also done by comparison with the exact SR-distance, 
that was computed by explicit formulas for SR-geodesics (given on p.386 in [25]) and the ex¬ 
plicit formulas for cut time (that coincides with the first Maxwell time, given by (5.18)-(5.19) 
in [25]). The experiments were done in the following way: 

1. Compute a set of points: 

EP(T) = {(xi,yi,0i) = Exp(pi,T) \ Pi E C,T < t^ AX (pi),i = 1,... i max } 

of end points lying on the SR-sphere of fixed radius T by analytic formulas for the 
exponential map and 1st Maxwell time t^ AX [25]. The number of end points i m ax 
was chosen as imax = 72 T 2 . Here C is the cylinder in momentum space given by 

c = {pe T:(SE( 2)) I Hf ixed (e,p ) = 1/2 } , 

where we recall (3.6). The sampling points pi are taken by a uniform grid on the 
rectifying coordinates (<£?, k) of the mathematical pendulum (the ODE that arises in 
the PMP procedure, cf. [25, ch:3.2]), both for the rotating pendulum case (pi G C 2 , 
yielding 5-curves) and the oscillating pendulum case (pi G Ci, yielding U- curves), 
where we note that C — C\ U C2. 

2. Evaluate the distance function = W°°(xi,yi,0i) obtained by our numerical 

PDE-approach in section 5 for every point of the set EP{T). We use 3rd order Hermite 
interpolation for W°°(xi,yi,Qi) at g = gi G EP(T ) in between the grid (5.1). 

3. Compute the maximum absolute error max I W(qA —T I and the maximum relative 

gi eEP(ry 

error max (\W (gA — T\/T). 
gi€EP(T) 

Remark 11. The exponential map Exp : C x M + SE( 2) provides the end-point 

(x(t),y(t),0(t)) = 7 (t) = Exp(po,t) of the SR-geodesic 7 , given SR-arclength t and initial 
momentum po G T*(SE( 2)). This exponential map integrates the PMP ODE’s in Appendix A. 

The absolute and relative errors of the SR-distance computations for each end points 
located on SR-spheres of fixed radii are presented in Figure 10. The red graph corresponds to 
a sampling of (TV, N$) = (50, 64), recall Eq. (5.1), used in the SR-distance computation by our 
numerical PDE approach, and the blue graph corresponds to the finer sampling (TV, Nq) = 
(140,128). We see that the maximum absolute error does not grow, and that the relative error 
decreases when increasing the radius of the SR-sphere. Increase of sampling rate improves the 
result. For the finer sampling case, neither the absolute errors nor the relative errors exceed 
0.1. 

6.1.3. Comparison and Computation of 1st Maxwell Set. Considering forward and 
backward derivatives acting on the distance function W°°(x,y,0), we can compute the 1st 
Maxwell set (recall eq. (3.5), see also Appendix D) as set of points where forward and back¬ 
ward left-invariant derivatives acting on distance function have different signs: 

2 

(6.2) M num = 1J {(x,y,0) € SE(2) \A+W°°(x, y, 9) > 0, ArW°°(x,y, 9) < 0}. 

2=1 
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Figure 9. Comparison of SR-spheres obtained by our numerical PDE-approach and the set of points 
EP(T ) lying on exact SR-spheres obtained by analytic formulas. From left to right: the SR-sphere with radius 
t — T — 3, T = 4 and T — 5. The color indicates the difference between the exact and the numerical values 
of the SR-distance (blue for smallest, green for middle, and red for highest differences). Thus, we see our 
algorithm is accurate, in particular along the fixed coordinate grid directions along x and 0-axis. 



A: Maximum absolute error B: Maximum relative error 


Figure 10. Maximum error in computing of SR-distance for end points located on SR-spheres of different 
radii t — T (from 1 to 7 with step 0.1), with number of end points imax = 72 T 2 . In red: errors are computed 
on a courser grid ( N , Nq) — (50, 64) and in blue: errors on a finer grid, with step sizes sq — and s xy — j-. 


Here i — 1 corresponds to the local component of the 1st Maxwell set (i.e. the part around e), 
and i — 2 corresponds to the global component of the 1st Maxwell set. The local component 
consists of two connected components lying on the surface given by xcos | + y sin | = 0 (i.e. 
the purple surface in Figure 5), and the global component is a plane given by equation 9 — n 
(for details, see [36]). In figure 11 we compare the local component of M nU m computed by 
our PDE approach with its exact counterpart M, presented in [36, thm 3.5]. It shows that 
Mnum is close to the exact Ai, except for the conjugate points at the boundary (where (6.2) 
is less accurate). Although note depicted here a similar picture was obtained for the global 
component, where M nU m indeed covers the plane 9 = tv. Summaring, this experiment verifies 
the correctness of the proposed method, but it also shows that the method allows to observe 
the behavior of the 1st Maxwell set. Eq. (6.2) allows us to numerically compute the Maxwell 
set for the data-driven cases C ^ 1 where exact solutions are out of reach. 

6.2. Feasibility Study for Application in Retinal Imaging. As a feasibility study for the 
application of our method in retinal images we tested the method on numerous image patches 
exhibiting both crossings, bifurcations, and low contrast, (Fig. 12, Fig. 13). For each seed 
point go the value function g W°° (g^ 1 g) was calculated according to the implementation 
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Figure 11. Comparison of the 1st Maxwell set obtained by our numerical PDE-approach with the exact 
1st Maxwell set [36]. Note that the local components of the 1st Maxwell set are part of the purple surface in 
Figure 5. Left: Local component of the exact Maxwell set M obtained by [36, thm 3.5] (where we recall that 
the cut locus coincides with the closure M of the first Maxwell set [35, th:3.3]). Middle: Local components of 
the Maxwell set M n um computed numerically by eq. (6.2). Right: Single case of a Maxwell point on the local 
part of the Maxwell set. 


details in Section 5, after which multiple end-points were traced back to the seed point. The 
image dimensions of the patches is 180 x 140. 

For the construction of the cost function (see e.g. Fig. 7) we set p = 3, and the lifting was 
done using cake wavelets with angular resolution tt/ 16. More precisely we used a cake-wavelet 
with standard parameters (N = 8, N$ = 32,5# = f ,cr s = 20px,7 = 0.8), for details see cf. [4, 
ch:2]. The precise choice of anisotropic wavelet is not decisive for the algorithm (so other type 
of anisotropic wavelets and cost constructions could have been applied as well). 

In all experiments we run with 4 settings for the two parameters (/3, A) determining the sub- 
Riemannian geodesics, we set p sma u = 0.05, fii arge = 0.1, X sm aii = 10, \arge = 100. The idea 
of these settings is to see the effect of the parameters, where we recall /? controls global stiffness 
of the curves, and A controls the influence of the external cost. We also include comparisons 
to a Riemannian wavefront propagation method on M 2 , and a Riemannian wavefront propa¬ 
gation method on SE( 2). These comparisons clearly show the advantage of including the sub- 
Riemannian geometry in the problem. For results on two representative patches, see Figure 12. 
For results on 25 other patches see www.bmia.bmt.tue.nl/people/RDuits/Bekkersexp.zip. 
Here, for fair and basic comparison of the geometries, we rely on the same cost function C. 
That is, we compare to 

• Riemannian geodesics y(t) = (x(t), y(£), 9{t)) in (S1E(2), G^ ull ) with 

C/J 7(t) (mm = (CMm 2 m)\ 2 +fi\m\ 2 +n 2 \m\ 2 ) 

• Riemannian geodesics x(s) = (x(s),y(s)) in (R 2 ,G^ 2 ) with metric tensor 

G m 2 Ix( s ) ( x ( s ), x ( s )) = (c(x(s))) 2 (|i;(s)| 2 + |y(s)| 2 ), 
with c(x(s),y(s )) = min ^C(x(s),y(s),0). 
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Typically, the wavefront propagation tracking methods on (M 2 , G^) produces incorrect short¬ 
cuts at crossings and very non-smooth curves. The Riemannian wavefront propagation track¬ 
ing method (with spatial isotropy) (SE(2),G ( j ull ) often deals correctly with crossings, but 
typically suffers from incorrect jumps towards nearly parallel neighboring vessels. Also it 
yields non-smooth curves. This can be corrected for when including extreme anisotropy, see 
Remark 12 below. The Sub-Riemannian wavefront propagation method produces smooth 
curves that appropriately deals with crossings. For high contrast images with reliable cost C 
best results are obtained with low (3 and large A. However, in low contrast images and/or 
patient data with severe abnormalities, low A is preferable, as in these cases the cost function 
is less reliable see Figure 13. 

Remark 12. It is possible to construct a family of anisotropic Riemannian metric tensors, 
recall (3.7)): G c e — C 2 (fPuo 1 0 go 1 + uj 2 0 uj 2 + /3 2 e 2 cu 3 0 uj 3 ), which bridges the SR-metric G c 
of our method (obtained by e ^ oo) to the full Riemannian metric tensor G C j ull (obtained by 
e 1). For the values of f3’s considered here, Riemannian geodesics and smooth Riemannian 
spheres for highly anisotropic cases e > 10 approximate SR-geodesics and non-smooth SR- 
spheres. In fact, with such extreme anisotropy in the Riemannian setting, the non-smooth 
ridges A4 in the SR spheres (see. e.g. the 1st Maxwell sets in Figure 4) are only little smoothed, 
and also the cusp-surface hardly changes. This observation allows to use the anisotropic fast¬ 
marching [24] as an alternative fast method for computing the solution of (3.1 ), instead of the 
iterative upwind finite difference approach in Section 4- 

The experiments indicate /? = 0.01 (small) in combination with A = 100 (large) are 
preferable on our patches. This typically holds for good quality retinal images of healthy 
volunteers. In lesser quality retinal images of diabetic patients, however, the cost function is 
less reliable and here A = 10 (small) can be preferable, see Figure 13. 

However, it might not be optimal to set the /? parameter globally, as we did in these 
experiments, as smaller vessels are often more tortuous and therefore require more flexibility, 
see e.g. [5, fig.7]. Furthermore, we do not include precise centerline extraction, which could 
e.g. be achieved by considering the vessel width as an extra feature (as in [7, 28, 21]). 

In future work we will pursue on sub-Riemannian fast-marching methods for fully auto¬ 
mated vascular tree extraction starting from an automatically detected optic nerve head via 
state-of-the art method [6] followed by piecewise sub-Riemannian geodesics (comparable to 
[14] in a Riemannian setting) in between boundary points detected by a ST?(2)-morphological 
approach. First experiments show that such a SR fast-marching method leads to considerable 
decrease computation time, hardly reduces the accuracy of the method, and can be used to 
perform accurate and fast automatic full retinal tree segmentation. The advantage of such an 
approach over our previous work on automated vascular tree detection [4] is that each curve 
is a global minimizer of a formal geometric control curve optimization problem. However, 
the SR-fast marching and automatic detection of the complete vascular tree via piecewise 
sub-Riemannian geodesics is beyond the scope of this theoretically oriented paper. 

7. Conclusion. In this paper we propose a novel, flexible and accurate numerical method 
for computing solutions to the optimal control problem (2.4), i.e. finding SR-geodesics in 
SE( 2) with non-uniform cost. The method consists of a wavefront propagation of geodesically 
equidistant surfaces computed via the viscosity solution of a HJB-system in (SE{2), A, G c ), 
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R 2 - Riemannian SE( 2) - Riemannian SE( 2) - Sub-Riemannian 



Figure 12. Data-adaptive sub-Riemannian geodesics obtained via our proposed tracking method 
(Thm. 3.2), with external cost (4-3), with p — 3, ft equals fismaii — 0.01, fiiarge = 0.1 and A equals 
A small = 10, A large = 100. We applied tracking from 2 seed-points each with several end-points (to test 
the crossings/bifurcations). To distinguish between tracks from the two seed-points we plotted tracts in different 
lighting-intensity. We indicated the valid cases only if all trajectories are correctly dealt with. 


and subsequent backwards integration, which gives the optimal tracks. We used PMP to 
derive both the HJB-equation and the backtracking. We have shown global optimality for 
the uniform case (C = 1) and that our method generates geodesically equidistant surfaces. 
Compared to previous works regarding SR-geodesics in (SE(2), A, G 1 ) [36, 18, 23], we solve 
the boundary value problem without shooting techniques, using a computational method that 
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Figure 13. Tractography (again for A = A small — 10, A = A large — 100 and — 0.01, p = 3) in a patch of 
a challenging low-contrast retinal image of a diabetic patient. In case of low-contrast (and less reliable cost) it 
is better to keep A small, in contrast to high contrast cases depicted in Figure 12. To distinguish between tracks 
from the two seed-points we plotted tracts in different lighting-intensity. We indicated the valid cases only if all 
trajectories are correctly dealt with. 


always provides the optimal solution. Compared with wavefront propagation methods on the 
extended domain of positions and orientations in image analysis [28, 29], we consider a SR 
metric instead of a Riemannian metric. Results in retinal vessel tracking are promising. 

Fast, efficient implementation using ordered upwind schemes (such as the anisotropic Fast 
Marching method presented in [24]) is planned as future work as well as adaptation to other 
Lie groups such as SE(3 ), SO( 3). Of particular interest in neuroimaging is application to high 
angular resolution diffusion imaging (HARDI) by considering the extension to SE( 3) [17, 29]. 
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Appendix A. Application of PMP for Canonical Equations for Cost-adaptive Sub- 
Riemannian Geodesics. We study optimal control problem (2.4). Recall Remark 2. Next 
we apply PMP to the action functional J Eq. (2.5) with fixed total time T > 0. Since 
[Ai,Aj] = Ylk=i c ij^k, with non-zero coefficients cf 2 = ~^ 2 i = ~ 1> c 23 = ~ c 32 = — 1, 
we have [A, A] = T(SE( 2)) and we only need to consider normal trajectories. Then the 
control dependent Hamiltonian of Pontryagin’s Maximum Principle (PMP) expressed via left- 
invariant Hamiltonians hi(p,g) = (p, Ai(g)), i = 1,2,3, with momentum p G T*(SE( 2)), and 
g — (x, y, 6 ) G M 2 x S 1 reads as 

H u (p,g) = v}-hi(p,g) +u 2 h 2 (p,g) - ^C 2 (g) (/3 2 \u l \ 2 + \u 2 \ 2 ) . 
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Optimization over all controls produces the (maximized) Hamiltonian 

and gives expression for extremal controls u 1 (t) = c2 ? u 2 (t) = ' Using SR- 

arclength parametrization G c (y(£), 7 (t)) = 1 implies H^ lxed = \ along extremal trajec¬ 
tories. We have the Poisson brackets 


(A.l) 


{H, M = 4r + {H, h 2 } = 


h i fe.s 
/3 2 C 2 ' 


{H,h 3 } = ^f 


/Z2^1 

C 2 ’ 


where H = W ixed and with {F, G} = £ $57 AG - fgA;F. By Eq. (A.l), by {hi,hj} = 

2=1 

3 

Aihj — Ajhi = X c ijhk, and by h t = {H, hi}, PMP gives us: 

k=l 


(A.2) 


(A-3) 


p( •) = J2 hi (' 


cj 


2=1 


l 7 (-) 


and < 


^ = A) Ai \a-) C+ c*$))’ 

= Al C ~ F^RO) ’ 

vertical part (for adjoint variables), 


* = c*mw cos6 ' 


7(1 = 53“*(-) A| 7 (.) and 2 y = sin6>, 

4 — ^2 
r _ c 2 (7(-))’ 

— horizontal part (for state variables). 


2=1 


with dual basis {cj 2 } for T*(SE(2)) defined by (uj l ,Aj) = 

For a consistency check, we also apply the PMP-technique directly to Problem (2.4) with 
free terminal time T, where typically [ 1 ] the Hamiltonian vanishes. Then, using SR arclength 
parameter t, the control dependent Hamiltonian of PMP equals 

H u (g,p) = v}hi{p,g) +u 2 h 2 (p,g ) - C(g)y /+ |u 2 | 2 . 


Optimization over all controls under SR arclength parametrization constraint 

l ^ 1 ! 2 + \u 2 \ 2 = 1 produces via EL-optimization w.r.t. (vf^u 2 ) (via unit Lagrange multi¬ 
plier) the (maximized) Hamiltonian: 

(A.4) H free (g,p) = + 1^21 2 - 1 = 0 with p = X 

and by straightforward computations one can verify that both the horizontal part and the 
vertical part of PMP (but now applied to H^ ree ) is exactly the same as (A3) and (A.2). 

Remark 13. The two approaches produce the same curves and equations, but their Hamilto¬ 
nians are different. Nevertheless, we have H^ ree = 0 O H^ lxed = 
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Appendix B. Lemmas Applied in the Proof of Theorem 3.1. 

In this section we consider preliminaries and lemmas needed for Thm. 3.1. Before we can 
make statements on SR-spheres we need to explain the notion of geodesically equidistant sur¬ 
faces, and their connection to HJB-equations. In fact, propagation of geodesically equidistant 
surfaces in (SE{2), A, G c ) is described by a HJB-system on this SR-manifold. 

Recall Remark 2 . Also recall, that in Appendix A we have applied PMP to this prob¬ 
lem yielding constant Hamiltonian H^ lxed — ^z{P~ 2 h\ + /i|) = \ relating to H^ ree — 
^^/Fy^Th% - 1 = 0 via Hf' ree = - 1. 

In our analysis of geodesically equidistant surface propagation we first resort to the non- 
homogenous viewpoint on the Lagrangian and Hamiltonian (with fixed time), and then obtain 
the results on the actual homogeneous problem (with free time) via a limiting procedure. 

Definition B.l. Given V : SE( 2 ) x M + —>► M continuous. Given a Lagrangian L(y(r), 7 (r)) 
on the SR manifold (SE(2 ), A, G c ), with L(y, •) : A —» M + convex. Then the family of surfaces 

(B.l) S r := {g G SE( 2 ) | V(g,r) = W 0 (r)}, with 

Wo :1gK monotonic, smooth, is geodesically equidistant if L(y(r), 7 (r)) = W((r). 

Remark 14. The motivation for this definition is 

R 

^ J L( 7 (r), 7 (r))dr = L( 7 (i 2 ), i{R)) = 

0 

Lemma B.2. Let L be non-homogeneous and lim = oo. Then the family of surf aces 

|v |—>oo ' V ' 

{5 r } rG M is geodesically equidistant if and only if V satisfies the HJB-equation (where r may 
be monotonically re-parameterized): 

(B. 2 ) W(g, r ) = -H(d SR V(g, r)), with d SR V(g , r) = F* A dV(g, r) = £ AV(g, r) . 

i —1 

3.2. 

Here P^(p = ^ hiuo 1 ) — ^ /q uo l is a dual projection expressed in dual basis d 1 given by 
2=1 2=1 

(d.Aj) — Sand Hamiltonian H('~v,p) = max {(p,v) — L(^y,v)}. 

J J ver^{SE(2)) lx 

Proof Substitute an arbitrary transversal minimizer y(r) into V(-,r) and take the total 

derivative w.r.t. r: 


f;V(j(r),r) = JA( 7 (r),r) + <dV| 7(r) , 7 (r)). 
Now 7(r) on S r , with tangent 7(r) = Ya =1 u% {t) and thereby 
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As a result we have 


(B.3) 


-i( 7 (r),i(r)) + £ u^r) A| 7 (r) E( 7 (r),r) = -f£( 7 (r),r) ^ 

i=1 

sup T,u i (r)h i (r) - T( 7 (r), 7 (r)) = -f£( 7 (r),r) Q 

(u 1 (r),u 2 (r))£l 2 i= 1 

H(F* A dV(j(r),r)) = 


with components hfir) = Ai\ 1 ^ V{^ (r),r) of projected 


momentum vector 


P( r ) W *l 7 (r) =I®A d ^(7(^),0 


2=1 


Now every point g G S r is part of a transversal minimizing curve y(r) and the result follows. 
So the “=4>” is proven. Conversely, if the HJB-equation is satisfied it follows by the same 
computations (in reverse order) that £( 7(^)5 7( r )) — (r),r), which equals W((r). □ 

Remark 15. In PMP [1] (see also Appendix A) the controls are optimized to obtain the 
Hamiltonian H from the control dependent Hamiltonian H u . The first equivalence in (B.3) 
is due to the maximum condition of PMP. The second equivalence in (B.3), is by definition 
of the Hamiltonian, where by the convexity assumption of the Lagrangian the supremum is 
actually a maximum [19, ch:8]. 

Next we apply the limiting procedure to obtain HJB-equations for geodesically equidistant 
surfaces in the actual homogeneous case of interest. The actual homogeneous Lagrangian case 
with T-free can be obtained as a limit (1 < g oo) from non-homogeneous Lagrangian cases: 


(B.4) 




(7 


2 rj V h(t) 

and corresponding Hamiltonian (see Remark 17 below) equals 


2ry —1 


(B.5) 


1 


H^{t),p{t)) = - {r 2 h\ + hiy ic( 7 (t))r 2 ", 


and setting r = t = W 0 (t). Thus f£( 7 (r),r) = %C({t),t) = W^t) = L( 7 (t), 7 (t)) = 
\jG C 1 7 ( t ) (j (2), 7 '( 2 )) = 1 in Eq. (B.3). Next we replace V by IT to distinguish between the 

homogeneous and the non-homogeneous case). 

Lemma B.3. The family of surfaces given by Eq. (B.l) is geodesically equidistant w.r.t. 

homogeneous Lagrangian £ 00 ( 7 , 7 ) = ^G c \ 1 ( 7 , 7 ), with r — t — Wo(t), iff W satisfies the 
HJB-equation 


(B. 6 ) 


ly/3- 2 |AhE | 2 + |^ 2 W | 2 = 1 ^ H = 0 


where H = lim H v = H^ ree the vanishing free-time Hamiltonian in Appendix A. Defining 

77 —>-oo 


Hamiltonian H by 
(B.7) 


H(i,p):=C 1 ( 7 ) JJ+hl 
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puts Eq. (B.6) in Eikonal form H(d SR W (g,t)) = 1 . 

Proof Tangential to the proof of Lemma B.2. For 1 < g < oo we can apply Lemma B.2 to 
Lagrangian L v given by (B.4) whose associated Hamiltion H] is given by (B.5) due to PMP 
(or just the Fenchel transform on M 2 ). In the limiting case g —>> 00 , where the Lagrangian is 
homogeneous and the Hamiltonian vanishes. Finally we note that now we have 

dW dW 

^(7 (r),r) = ^-(7 = W'(t ) = L( 7 (t), 7 (t)) = 1 , 

from which the result follows. □ 

Remark 16. The relation between the various Hamiltonians is 

oo = H free = - 1 = y/2H v= i -1 = H-1 = 0. 


Remark 17. The relation between the Lagrangian L v given by (B.4) and the Hamiltonian 
(B.5) is the (left-invariant, SR) Fenchel transform on SE( 2). Due to left-invariance this 
Fenchel-transform actually boils down to an ordinary Fenchel-transform on M 2 when expressing 
velocity and momentum in the left-invariant frame. Indeed we have 

(B.8) H v (^,p) = [Sc(SE(2))nA(L r] (j, ■))](?) : = 

sup {———-(C(7)) ^FT (/ 3 2 |t/ 1 | 2 + |iff | 2 ) 277—1 _|_ h lU l h2U 2 } 

(u 1 ^ 2 )^ 2 2rf 

with horizontal velocity v = u l A\ + u 2 A 2 , and momentum p = Y)n=i hiU 1 . 

Appendix C. Viscosity Solutions for HJB-systems in SE(2). 

Definition C.l .The (Cauchy problem) for a HJB-equation (akin to [19, ch:10.1]) on SE( 2) 
is given by 

f^ = -H(g,d SR W) in SE( 2) x (0, T), 

\W(g,0) = Wo, 

whereas a boundary value problem for HJB-equation is given as 
(C.2) H(g ,= 0 on SE( 2) \ {e}, W(e) = 0; 

where T > 0 is prescribed, Wq is a given function (or a cost measure [2]), H(g,p) = H-f ree (g. p) 

2 

is the free-time Hamiltonian given by (3.8), and d SR W = ^ AiW(g,t) cv l \ . 

Remark 18. Combined Cauchy-Dirichlet problems exist [38], but they are defined on (an¬ 
alytic) open and bounded domains. Thereby they cannot be applied to our set of interest 
SE{ 2) \ {e} as this would violate semigroup theory [2, 19, 34, 12]. This is also clear in view 
of the Cramer transform [2], putting an isomorphism between HJB- and diffusion systems. 

Remark 19. In Eq. ( C.2) it is crucial that the free time Hamiltonian is used. In the definition 
of viscosity solutions of the Cauchy problem, Eq. (C.l), one can set both H = H^ ree (as done 
in the body of the article) or H — H^ lxed as done in Appendix B. 
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HJB-systems in general do not have unique (weak) solutions. To avoid multiple (non¬ 
desirable) solutions, one must impose the viscosity condition [19, 22] commonly applied in 
wavefront methods acting directly in the image domain M 2 [27, 37]. The viscosity solution is 
obtained by the vanishing viscosity method [22]. The idea of this method is to add to the 
HJB-equation a term e A and to pass to the limit, when s goes to 0. Here A denotes the 
Laplacian, that in our case (for C — 1) equals A SR = Y^i=i A(A) -2 A? with /?i =/?,/? 2 = 1- 
Here the name “viscosity solutions” comes from fluid dynamics, where typically the term 5 A 
represents a physical viscosity. For an intuitive illustration of the geometric property of such 
solutions see [11, fig.30]. The viscosity solution of the initial value problem can be defined 
alternatively as follows. 


Definition C.2 .Let H(g , •) be a convex Hamiltonian for all g G SE( 2) s.t. H(g,p) —)> 00 if 
p —>► 00. The function W : SE( 2 ) x M —>► M is viscosity solution of = —H(g,d SR W) if it 
is a weak solution such that for all functions V G C 1 (SE( 2) x M,M) one has 

• ifW — V attains a local maximum at (go, to) then ®+j/( 9 ,d^))i < 0 , 

• ifW — V attains a local minimum at (go, to) then (^ + ff( 9 ,d^))| >0. 

Similarly, the viscosity solution of the boundary value problem (that is equivalent to 

Eikonal equation, when t is SR-arclength) can be defined as follows: 

Definition C.3. A solution W : SE( 2) -G M of Eq. (C.2) is called a viscosity solution if for 
all functions V G C 1 (SE{2),M) one has 

• ifW — V attains a local maximum at go then Hf ree (go,d SR V)) < 0 ; 

• ifW — V attains a local minimum at go then Hf ree (go, d^R)) > 0. 

Appendix D. Proof of Theorem 3.2 . 

The back-tracking (3.3) is a direct result of Lemma B.3 in Appendix B and PMP in Ap¬ 
pendix A. According to these results one must set 


u\t) = 


hi(t) 

( C ( 7 ( t ))) 2 /? 2 


EEEr an( j u 2(t) = = EEEri 

(c( 7 (i))) 2 /3 2 1J {cm )) 2 {cm )) 2 


from which the result follows. Then we recall from Thm. 3.1 that St given by (3.2) are 
geodesically equidistant surfaces propagating with unit speed from the origin. So St are 
candidates for sub-Riemannian spheres, but it remains to be shown that the back-tracking 
(3.3) will neither pass a Maxwell point or a conjugate point, i.e. t < t cu t. Here t cu t denotes 
cut time, where a geodesic looses its optimality. 

At Maxwell points g * induced by the 8 reflectional symmetries [25] two distinct SR- 
geodesics meet with the same SR distance. As SR-geodesics in (SE(2), A, G 1 ) are analytic 
[25], these two SR-geodesics do not coincide on an open set around end-condition g*, and the 
SR spheres are non-smooth at g*. Regarding the set Ai, we note that the Maxwell sets related 
to the i-th reflectional symmetry C{ are defined by 

MAX* = {(po,£) e T*(SE( 2)) x M+ | H(p 0 ) = \ and Exp{p 0 ,t) = Exp{eip 0 ,t)}, 
max z = Exp( MAX Z ), z = 1,... 8 , 


where we may discard indices i = 3,4, 6 as they are contained in {maxi, max 2 , max 5 , max 7 }. 
Now with max 2 we denote the Maxwell set with minimal positive Maxwell times over all 
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symmetries (i.e. we include the constraint t < min {t l max } where the minimum is taken over 
all positive Maxwell times along each trajectory), then we find A4 to be contained within the 
union of the following sets 3 * : 

max 2 C {(x,y,0) G SE( 2) | y sin 9/2 + x cos 9/2 = 0}, max 5 = {(x,y,0) G SE{2 ) | 9 = i r}, 

where [25, th:5.2] shows we must discard the first reflectional symmetry e\ as it does not 
produce Maxwell points. Now for generic geodesics (not passing the special conjugate points 
that are limit points of Maxwell points and not Maxwell points themselves) t cu t — as 

proven in [35, th:3.3], where t x MAX > 0 denotes the first Maxwell time associated to the 8 
discrete reflectional symmetries. 

During the back-tracking the set Ai is never reached at internal times (only when starting 
at them, recall Remark 6), since they are “uphill” from all possible directions during dual 
steepest descent tracking (3.3), as we will shown next. As a result we have t < t cu t — t x MAX . 
Consider Fig. 5. At Maxwell points y* G A4 due to the reflectional symmetries there exist two 
distinct directions in the 2D-horizontal part A g * of the tangent space T g *(SE( 2)) where the 
directional derivative is positive. If there would be a direction in the tangent space where the 
directional derivative is negative then there would be a direction in A g * with zero directional 
derivative of W(-) at y* towards the interior of the sphere yielding contradiction. Here we 
note that due to the viscosity property of the HJB-solution, kinks at the Maxwell points are 
pointing upward (see Fig. 5 and Fig. 14) in the backward minimization tracking process [11, 
fig.30]. Furthermore, we note that SR spheres St are continuous [36] and compact, as they 
are the preimage St = d(-, e)^“({t}) of compact set {t} under continuous mapping d(-,e). 
Continuity of d(-, e) implies the spheres are equal to the 2D-boundaries of the SR balls (w.r.t. 
the normal product topology on R 2 x S' 1 ). 

The algorithm also cannot pass conjugate points that are limits of 1st Maxwell points, 
but not Maxwell points themselves. See Fig. 4. Such points exist on the surface R 2 — 0 and 
are by definition within M\M. Suppose the algorithm would pass such a point at a time 
t > 0 (e.g. there exist 4 such points on the sphere with radius 4, see Fig. 5) then due to the 
astroidal shape of the wavefront at such a point, cf. [35, Fig.11], there is a close neighboring 
tract that would pass a 1st Maxwell point which was already shown to be impossible (due to 
the upward kink-property of viscosity solutions). □ 

Remark 20. The sub-Riemmannian spheres are non-smooth only at the 1st Maxwell set A4. 
They are smooth at the conjugate points in M\M (where the reflectional symmetry no longer 
produce two curves/fronts). In the other points on St \M the SR-spheres are locally smooth 
(by the Cauchy-Kovalevskaya theorem and the semigroup property of the HJB-equations). 

Appendix E. The Limiting Procedure (4.3) for the sub-Riemannian Eikonal Equation. 

In this section we study the limit procedure (4.3), illustrated in Figure 6. To this end we 
first provide a formal representation of the viscosity solutions of system (4.1), where we rely on 
viscosity solutions of morphological scale spaces obtained by super-position over the (min, +) 

3 In [36, Eq.3.13] it is shown that max 2 = {(x,y,6) G SE( 2) | y sin 6/2 + x cos 6/2 — 0 and | — xsm{0/2) + 

ycos(#/2)| > |R}(0)|} with R\ defined in [36, Lemma 2.5]. We also observed such a loss of the Maxwell 

point property in our numerical algorithm, as kinks in W(g) = t can disappear when moving on the set 
y sin 0/2 + x cos 6/2 — 0. See Figure 11 . 
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Figure 14. Overview of Maxwell points. Two Maxwell points on the purple surface xcos | + y sin | = 0 
and two red Maxwell points on the surface \6\ = tt (recall Figure 5). In all cases we see that local kinks in the 
viscosity solutions are upward, and the back-tracking algorithm can not pass these points. 


algebra, i.e. obtained by morphological convolution (erosions) with the morphological impulse 
response, cf. [12]. Now as the HJB-equations of such morphological scale spaces do not involve 
a global offset by 1 in the right-hand side of the PDE we need to combine such erosions with 
a time shift in order to take into account the global offset. It turns out that the combination 
of these techniques provide staircases with steps of size e, so that we obtain the appropriate 
limit by taking the limit e —>> 0 afterwards as is done in (4.3). 

Morphological convolutions over the SE( 2) group are obtained by replacing in linear left- 
invariant convolutions (likewise the SE(3)-ca,se studied in [17]), the usual (+, *)-algebra by 
the (min,+)-algebra. Such erosions on SE( 2) are given by 

(E.l) (k © f)(g) := inf {k{h~ l g) + f(h)}. 


Furthermore, to include the updating of the initial condition in (4.1) we define 

< E2 ) « , W : ={r <9) if/=e e 


Lemma E.l. Let e > 0, n G N. The viscosity solution of (4-1) is given by 
(E.3) WZ +1 (g, r) = ( k r . ne © W$(g) + (r - ne), 

for r G [r n ,r n+ i] = [ne, (n + l)e] and morphological kernel k v (g), v > 0, given by 

h(g) = l° 

w oo else. 


where d(g,e) denotes the Carnot-Caratheodory distance (3.4) between g G SE( 2) and e = 
(0,0,0). Forn — 0 we have that the viscosity solution of (4-3) is given by W((g,r) — k r (g)+r. 

Proof In order to care of the constant off-set in the HJB-equatons of (4.1) and (4.2), we 
set r = r new + r n and we define for n = 0,1, 2, ... the functions V^ +1 : SE( 2) x [0, e] —>> R by 


Vn+1 (g,r n 


:= W^ +1 (g, r n 


4 “ r n ) r n 


with r new G [0, e] and V^ +1 the viscosity solution of 

dr)tw (g 5 rnew ) = — 1 4 " 1 — H(d SR Vn+i(g,r new )) — —H{d SR V^ 1 {g^r new )) J 

< for g ^ e we have V* + Ag, 0) = i \ !~ U f’ 

| W*(g,r n ) if n € N 

for g = e we have V* +1 (g, 0) = V^ +1 (e, 0) = 0. 
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where we use short notation for the sub-Riemannian derivative d SR V := reca ll 

(3.7) in Remark 5, and where Hamiltonian H is given by (B.7). 

Now let us first consider the case n — 0. By the results in Appendix B the Hamiltonian 
system (4.2) provides geodesically equidistant wavefront propagation traveling with unit speed 
and departing directly from the unity element.As a result, we find 

yet ) — k (n) — I ® ^ d(g,e ) < r new 

Vi (, 9 , rnew) - k rnew (g) - j ^ elge 

and by left-invariant ‘superposition’ over the (min,+)-algebra we find for n — 1,2,... that 
V n+\(g,r n ew) = ( K new © W*(-,r n ))(g). Finally, we have 

w n+i(9,r) = Vf +1 (g,r -ne) + r - ne = (fc r _ n€ eW„(-,r n ))(g) + r - ne. 


Corollary E.2. Let n G N. Let e > 0. The following identity holds 
W^ +1 (g,r n+1 ) = (k e eWfr,r n ))(g)+e 
(E-4) f X (m + l)e 1 [r m ,r m+1 ](d(g, e)) if d(g, e) < r n+1 = (n + l)e 

\ m =0 

[ oo if d(g,e) > r n+ 1 , 

where l[r m ,r m+ i] denotes the indicator function on set [r m ,r m +i]. 

Proof The first part follows by Lemma E.l for r = r n+ 1 (i.e. r new = e). The second part 
follows by induction. Recall from Lemma E.l that W[(g, r) — k r (g) + r. Now application of 
(E.3) for n — 1 yields 

Wi(g,r 2 ) = (k e eWi(;r 1 ))(g) + e = e +^^ k f h) + e * Jj ^ = 

+1 ]{d{g,e)) if d(g, e) < r 2 
else, 


(E.5) 


0 if d(g, e) < e 
e if e < d(< 7 , e) < 2e = 
oo else 


X ( TO - 

m— 0 

oc 


l)el [r 


with Bg :C = {h G SE( 2) | d(g,h) < e}. This can intuitively be seen from the geometric 
meaning of an erosion VEf k e © VEf where one drops cylinders from below on the graph of 
(-,r n ) and considering the new hull where cilinders get stuck. Eq. (E.5) can also be seen 
directly from the definition of k e . Let us verify each case separately: 

• If d(g,e) > 2e we have that the value must be infinite, since suppose it were finite 
then by the definition of the morphological kernel k e we would need to have that 
d(g , e) < d(g , h ) + d(h , e) < 2e yielding contradiction. 

• If d(g, e) < e, then in the erosion-minimization we can take h — e and we obtain e + 0. 

• If e < d(g , e) < 2e, then in the erosion-minimization we cannot take h = e, but for 
allowed choices we obtain k e (e) — 0 and e + e as output. 

Similarly we have by inserting induction hypothesis for n and recursion (E.3) we have 


Wf +2 (.g,r n+2 ) = {k e QW* +1 (-,r n+1 ))(g)+e 


Tl-|-1 

e+ X (™ + !) e ![r m+ 1 j7"m+ 2 }(d(g,e)) 

171=0 


TL ~h 2 

= e+ X ^ ,el [r- m „r m , +1 ](%,e)) 

m’=l 


n+1 

X K + l)el[r m „r m , + 1 ](%,e)), 

m'=0 
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for d(g , e) < r n + 2 . Here we applied m! = m + 1 so that the result follows for n + 1. 
Theorem E.3.Le£ g E SE(2) be given. We have the following limit 


lim lim W* +1 (0, (n + l)e) = d{g, e ). 

e —)>0 n—>• 00 


□ 


Proof Application of Corollary E.2 gives 

00 JV*(«/,e) 

lim W' +1 (fli, (n + l)e) = E ( k + U l [r k ,r k+1 ] (<% 0) = E (k + l)e l [r r ](d{g, e)), 

n ^°° k =0 k =0 

with N*(g,e ) = i.e. the largest integer > € M + . As a result we have 

N*(g,e) 

lim lim W* +1 (g,r n+1 ) = } +1 (g, (n + l)e) = lim E) ( fc + l) e \r k ,r k+l ]{d{g, e)) = d{g, e) 

where the size of the steps in the staircase towards d(g , e) vanishes as e —» 0. Recall Figure 6. 
Together with Theorem 3.2 yielding W°°(g) = d(^,e), the result follows. □ 
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